%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time-dependent optical Bloch equations for a two-level system,
%  n = [11 22 12 21],
%
%  --- |2> 
%	| 
%  --- |1> 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Bloch
    figure(1);			clf;
    col='rgbmc';

    G12=   1;                                                               % decay of population
    O12=[0:.05:1]*G12*8;                                                    % Rabi frequency
    d12=   0;                                                               % detuning

    n0 =[1 0 0 0]';                                                         % initial condition
    Tim=[0:.05:6];
    for kk=1:length(O12)
        o12=O12(kk);
        for tt=1:length(Tim)
            t = Tim(tt);
            n = Bloch2(n0,t,G12,o12,d12);
            nt(kk,tt)=n(2);
        end
    end
    subplot(221),mesh(O12/G12,G12*Tim,nt');
    set(gca,'FontSize',16)
    xlabel('\fontsize{16}\Omega_{12}/\Gamma_{12}');
    ylabel('\fontsize{16}\Gamma_{12}{\itt}');
    zlabel('\fontsize{16}{\it\rho_{22}}');
    text(1,-1.5,1.4,'\fontsize{18}(a)')
    view([75 40])

% sequence

    dt= .1;                                                                 % time step
    G = .1;                                                                 % decay rate
    O =  2;                                                                 % Rabi frequency
    D = [0 1];                                                              % detuning
    T = pi/O;                                                               % pi-pulse duration
    
    for kk=1:2
       n =[1 0 0 0]';
       d =D(kk);
       tt=    0;
       for tt=0:dt:T+20
            if tt<T
                o = O;          g = G;                                      % pi-pulse
            else
                o = 0;          g = G;                                      % spontaneous decay
            end

            ns(1)=sqrt(2)*imag(n(3));
            ns(2)=sqrt(2)*real(n(3));
            ns(3)=real(n(2)-n(1));
            n = Bloch2(n,dt,g,o,d);
            nt(1)=sqrt(2)*imag(n(3));
            nt(2)=sqrt(2)*real(n(3));
            nt(3)=real(n(2)-n(1));

            subplot(222),plot3([ns(1) nt(1)],[ns(2) nt(2)],[ns(3) nt(3)],col(kk),'MarkerSize',10,'LineWidth',1.2)
            set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',16,'LineWidth',1.2)
            xlabel('\fontsize{16}Im {\it\rho}_{12}')
            ylabel('\fontsize{16}Re {\it\rho}_{12}')
            zlabel('\fontsize{16}2{\it\rho}_{22}-1')
            text(-1,-1.5,1.5,'\fontsize{18}(b)')
            view([75 40])
            hold on
            pause(.01);
        end
    end
    grid
    hold off
end

function psi = Schrd2(psi0,t,o,d)
    H =[0    o/2 
        o/2   d ];
    psi=expm(1i*H*t)*psi0;
end

function n = Bloch2(n0,t,g,o,d)
    I =     1i/2;

    M=[ 0      g     I*o   -I*o 
        0     -g	-I*o    I*o 
        I*o  -I*o -g/2-1i*d  0   
       -I*o   I*o     0   -g/2+1i*d];
    n=expm(M*t)*n0;
end
